Entanglement Spectra of the 2D AKLT Model: VBS/CFT Correspondence 
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We investigate the entanglement properties of the valence- bond- solid (VBS) state defined on two- 
dimensional lattices, which is the exact ground state of the Affleck-Kennedy-Lieb-Tasaki model. 
It is shown that the entanglement entropy obeys an area law and the non-universal prefactor of 
the leading term is strictly less than In 2. The analysis of entanglement spectra for various lattices 
reveals that the reduced density matrix associated with the VBS state is closely related to a thermal 
density matrix of a holographic spin chain, whose spectrum is reminiscent of that of the spin-1/2 
Heisenberg chain. This correspondence is further supported by comparing the entanglement entropy 
in the holographic spin chain with conformal field theory predictions. 



PACS numbers: 75.10.Jm, 75.10.Kt, 05.30.-d, 03.67.Mn 



I. INTRODUCTION 

There has been considerable recent interest in 
understanding quantum entanglement in many-body 
systems*^ Entanglement (von Neumann) entropy and 
the family of Renyi entropies are often used to quantify 
the degree of entanglement of a bipartite system consist- 
ing of two subsystems A and B. However, as pointed out 
by Li and Haldane^ 3 - the entanglement spectrum (ES), 
the eigenvalue spectrum of the reduced density matrix 
(RDM) of either system A or B, provides a more com- 
plete description of the ground state and the low-energy 
states of the system. In general, the RDM for A may be 
written as pa = exp(— 7?e) with the entanglement Hamil- 
tonian He- Li and Haldane demonstrated that the spec- 
trum of -He for fractional quantum Hall states on a sphere 
reflects the gapless edge excitations in the disk geometry. 
Since then, the ES has been applied to other systems 
such as topologically ordered systems^— and quantum 
spin models^— and the entanglement gap was proposed 
as a non-local order in gapless spin chains^ 

Although the analysis of the ES has turned out to 
be a powerful tool to extract universal properties, the 
application has so far been limited to one-dimensional 
(ID) or topological systems. To fill this gap in the lit- 
erature, we examine the ES in the two-dimensional (2D) 
quantum spin model proposed by Affleck, Kennedy, Lieb, 
and Tasaki (AKLT). -±2. The AKLT model in one dimen- 
sion was first introduced as a solvable model to illustrate 
the Haldane-gap phase in integer spin chains. The exact 
ground state is known as the valence-bond-solid (VBS) 
state. The construction of the VBS state generalizes to 
two and higher dimensions i 14 ^ 15 The VBS state has re- 
cently attracted renewed interest from the viewpoint of 
quantum computation. It was proposed that VBS states 
on a 2D hexagonal lattice can serve as resources for 
measurement-based quantum computation i 16 i Thus, 
the study of the ES in VBS states on various 2D lattices 
can shed some light on this important issue. 

In this paper, we study numerically the ES in the 2D 



AKLT model defined on various lattices. We develop an 
efficient way to evaluate the RDM using a Monte Carlo 
(MC) method. The ES is then obtained by numerical 
diagonalization. It is shown that the ES is not flat in 
2D, which is in contrast to the ID case where all eigen- 
values of the RDM become degenerate in the thermo- 
dynamic limit i 18 i 19 When a system is on a cylinder, we 
can study the ES as a function of momentum. The low- 
lying spectrum of the entanglement Hamiltonian He for 
square (hexagonal) VBS state is reminiscent of that of the 
ID antiferromagnetic (ferromagnetic) Heisenberg chain. 
This indicates that the entanglement Hamiltonian He is 
very close to the Heisenberg Hamiltonian in one dimen- 
sion. To further support this idea, we introduce a new 
concept, nested entanglement entropy (NEE), and find 
that, for the square lattice VBS state, the low-energy 
effective Hamiltonian of He is well described by c = 1 
conformal field theory (CFT) despite that the starting 
quantum state (the AKLT state in this paper) is not 
critical. For a hybrid system consisting of squares and 
hexagons, we find that the ES is gapped and He is close 
to the alternating Heisenberg chain in which ferromag- 
netic and antiferromagnetic interactions are mixed. 

The organization of the rest of the paper is as follows. 
In Sec. II, we first review the construction of the VBS 
ground states on two-dimensional lattices. We then show 
how to obtain the spectrum of RDM from the Schmidt 
decomposition of the VBS state. We see that the Schmidt 
coefficients are related to the eigenvalues of the overlap 
matrix (see Eq. ([6]) for the precise definition). In Sec. Ill, 
we show numerical results of entanglement entropy and 
spectrum of the VBS states on both square and hexag- 
onal lattices. We then introduce the NEE and study its 
scaling properties. We also present an example of the 
gapped entanglement spectrum, which is obtained from 
the VBS state on the hybrid lattice consisting of squares 
and hexagons. We conclude with a summary in Sec. IV. 
Analytical and numerical approaches for obtaining the 
overlap matrix are, respectively, presented in Appendices 
A and B. 



II. 2D VBS STATES 



B. Reduced density matrix 



A. AKLT Hamiltonian 

We study the AKLT model on the square and hexago- 
nal lattices shown in Fig. Q] (a). For simplicity, we focus 
on the basic model in which there is a single valence bond 
on each edge. The spin-S* Hamiltonian for this model is 
written as a sum of projection operators as 



H = ^2A{i,j)TT 2S {i,j), 



(1) 



(id) 



where A(i,j) is an arbitrary positive number and 
7i"2s(J,i) projects the total spin of each pair of neigh- 
boring spins onto the subspace of J = 25. The spin 
magnitude 5 and the coordination number z are related 
through z — 25, where z = 4 for the square lattice, while 
z = 3 for the hexagonal lattice. The unique ground state 
(VBS state) of H is written as 



\^) = \{{a\b]-b\a])\^c) 

(id) 



(2) 



where a] and b\ are Schwinger bosons at site i, satisfying 



[ai,a\] = [bi,b]} = S it . 



(3) 



with all other commutators vanishing, and |vac) denotes 
the vacuum stated With the constraint that the total 
number of bosons at each site is 25, the original spin 
operators can be represented by the Schwinger bosons as 



St = dbi, S: 



bleu, 5? 



-(a\ai-b\h). 



(4) 



We will study VBS states defined on lattices with both 
periodic (PBC; cylindrical geometry) and open (OBC; 
rectangular geometry) boundary conditions^ shown in 
Fig. [1(b). 



subsystem A , subsystem B 
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Let us consider a bipartition into two subsystems A 
and B, each of which is a mirror image of the other, and 
the boundary between them is chosen to be the reflection 
axis. We denote the bulk width of A (B) by L x and the 
number of bonds on the boundary by L y . The Hamilto- 
nian can be written in the form H — Ha + Hb + Hab , 
where Ha and Hb denote the Hamiltonians within sub- 
systems A and B, respectively, while Hab contains all 
the interactions across the boundary. The RDM for A 
describing the entanglement between the two subsystems 
is defined by p A = Tr B [|*)(*|]/A/', where Af = (W|W). 
To obtain the spectrum of pa, we follow the same ap- 
proach used in Ref. (22|. The ground state of the system 
can be written as 



i*) = Ei^ ] ) 



(5) 



where a = 1,2,...,2 L » and -j>L A] )} ({|4? ] }}) is a set of 
degenerate ground states of Ha (Hb), each of which is 
characterized by a configuration of bosons at the bound- 
ary sites. These states are linearly independent but not 
orthogonal. We now introduce overlap matrices as 



Ml A J = (<^W) 



l a/3 



M 



[B] _ 

a/3 ~ 



(6) 



for two subsystems, respectively. Due to the reflection 
symmetry, one can show that M^ = M ^ = M and M 
is real symmetric and positive definite. One can then con- 
struct the orthonormal states out of {|^L )} {{\<t>a )}) 
using an orthogonal matrix O that diagonalizes M, i.e., 
O t MO — diag(di,d 2 , ...,d 2 L y ): 



[A]\ 



\rl>. 



14 s 



<a a p 



(7) 



The ground state of the full system in the new basis is 
written in the Schmidt decomposition as 



l*) = E d -l^ 1 )®!^ 



(8) 



and hence pa — J^a^al^a }(wa |/A/". Therefore, the 
spectrum of pa is equivalent to that of the following ma- 
trix 



FIG. 1: (Color online) (a) VBS states on the square and 
hexagonal lattices with reflection symmetry. Lattice is sep- 
arated into two subsystems A and B by the reflection axis 
(the dotted line). The bulk width is denoted by L x , while the 
number of bonds on the boundary is L y . Notice that in the 
hexagonal lattice, L x is defined as half the number of sites 
along the x-axis. In both cases, L x = 3 and L y — 4. (b) 
Boundary conditions. Upper and lower panels show cylindri- 
cal (PBC) and rectangular(OBC) geometries, respectively. 



M 2 



PA 



Tr[M 2 ] ' 



(9) 



where Tr denotes the standard matrix trace (not to be 
confused with Tr^). One can regard pa as a thermal 
density matrix of an auxiliary spin chain, which we call 
a holographic spin chain. The notion of this hidden spin 
chain was first introduced in Ref. 1221 The entanglement 
Hamiltonian is defined via pa — exp(— He)- 



We can obtain the overlap matrix M analytically for 
VBS states defined on vertical ladders (i.e. L x — 1) (see 
Appendix |A|) . For VBS states on lattices with a larger 
width L x (> 1), in which case results converge to the 2D 
limit, we develop an efficient algorithm of MC method 
(see Appendix [B] for details). We sample all possible 
configurations of bosons stochastically, and each matrix 
element M a p can be numerically obtained as the accu- 
mulation number of states. 



III. ENTANGLEMENT ENTROPY, 

SPECTRUM, AND NESTED ENTANGLEMENT 

ENTROPY 
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In this section, we show the results of entanglement en- 
tropy and spectrum for the VBS states obtained by the 
MC method combined with exact diagonalization. As we 
will show, the ES of the square (hexagonal) VBS state 
resembles the energy spectrum of the antiferromagnetic 
(ferromagnetic) Heisenberg Hamiltonian in one dimen- 
sion. We also introduce a new quantity, nested entan- 
glement entropy, and further elucidate the relationship 
between the holographic spin chain and the Hamiltonian 
for the Heisenberg chain. 



A. Entanglement entropy 

Let us first study the entanglement entropy (EE) which 
can be obtained from 



S = -Tr[pA In pa] = ~^Pa lnp Q , 



(10) 



where p a (a = 1, 2, ..., 2 Ly ) are the eigenvalues of pa- 
The obtained EE divided by the boundary length for the 
square and hexagonal lattices are shown in Fig. [2] In 
both cases, we investigate lattices with both PBC and 
OBC, results of which are expected to be equivalent in 
the thermodynamic limit (L v — > oo). The EE per unit 
boundary length S/L y in 2D AKLT model is found to be 
strictly less than In 2 = 0.693147, which is in contrast to 
the result for the ID AKLT modeli±&i2i22, Note that the 
area law is satisfied, as S/L y approaches to a constant in 
the thermodynamic limit. 

The difference between ID and 2D EE results can be 
understood in terms of "valence bond loops" (VBL), 
which are formed by sequences of singlets connected 
head-to-tail. (A precise definition of VBL in the con- 
text of Schwinger boson can be found in Appendix FBI ) 
The existence of such loops is a unique feature of 2D 
VBS states, which provides additional correlations be- 
tween two separated singlets. As a result, we expect 
extra correlations between boundary bosons in the 2D 
VBS states, which lead to modification of the EE from 
the ID result. The possibility of forming a valence bond 
loop decays exponentially fast with respect to its length. 
(The reason is discussed in Appendix [Bj. This can be 



FIG. 2: (Color online) Entanglement entropy 5 per valence 
bond across the boundary as a function of L y for (a) square 
lattices with rectangular geometry (OBC), (b) square lattices 
with cylindrical geometry (PBC), (c) hexagonal lattices with 
OBC, and (d) hexagonal lattices with PBC. A further increase 
of L x does not affect the results much. The extrapolations to 
infinite size are indicated by the broken lines. 



related to multiple observed pheonomena of EE from our 
numerical calculation. 

First of all, in VBS states on the square lattice, the 
magnitude of deviation from the ID result {S/L y = In 2) 
is significantly larger than that of the hexagonal lat- 
tice. To understand this, we observe that elemental loops 
(shortest VBLs that have largest possibility) connects 
two boundary sites directly in the square lattice, while an 
extra site is in between two boundary spins in the hexag- 
onal lattice. As a result, the amplitude of elemental loops 
is higher in the square lattice. In the end, we obtain con- 
siderably more VBLs, which then induce larger amount 
of extra correlations on the boundary in the square lattice 
than in the hexagonal lattice. 

Meanwhile, increasing the number of sites L x in the 
bulk direction provides additional valence bond loops 
which involve sites not only on the boundary, but also 
inside the two-dimensional lattice. Hence the EE is fur- 
ther reduced when expanding bulk width L x . The effect 
converges exponentially fast, however. The 2D limit has 
been almost reached at L x — 5 and L x — 3 for square 
and hexagonal VBS states, respectively. It is again re- 
lated to the fact that the amplitude of a VBL decreases 
with its length. 

We also observe that the thermodynamic limit (L y —> 
oo) of EE is achieved with totally different behaviors for 
lattices with rectangular (OBC) and cylindrical (PBC) 
geometries, as shown in Fig. [2J For VBS states on a lat- 
tice with OBC, S/Ly is found to scale with the boundary 
length {L y ) as 



S_ 
1~, 



Co 



C x 



Ly Ly fog Ly 



(11) 



where C\,C% (> 0) are the fitting parameters and a 
denotes the EE per boundary length in the thermody- 



namic limit. The magnitude of Ci is significantly less 
than Ci, as shown in Table [J] In contrast, for a lattice 
wrapped on a cylinder (PBC), the thermodynamic limit 
is approached exponentially fast in L yi i.e., a — S / L y oc 
exp(— L y /£), where £ (also shown in the table) is of the 
order of the lattice spacing. 





square 






OBC 


PBC 




a C\ Ci 


t 


L x = 2 
L> x — o 

L X =A 
L x = 5 


0.61277(3) 0.0865(3) 0.0009(2) 
0.59601(4) 0.1127(4) 0.0034(3) 
0.59322(5) 0.1196(5) 0.0049(3) 
0.59259(5) 0.1223(5) 0.006(3) 
0.59246(4) 0.1227(5) 0.006(3) 


0.6129(1) 1.54(5) 
0.5966(2) 2.4(1) 
0.5942(3) 2.9(2) 
0.5936(4) 3.1(2) 
0.5934(4) 3.1(2) 





hexagonal 






OBC 


PBC 




<7 Ci C2 


t 


L« = l 

L x = 2 
L x = 3 


0.68502(2) 0.0092(6) 0.0011(7) 
0.68469(3) 0.0098(6) 0.0012(8) 
0.68468(3) 0.0098(6) 0.0013(7) 


0.68508(6) 0.4(1) 
0.684757(4) 0.6(1) 
0.684754(5) 0.7(1) 



TABLE I: Obtained fitting parameters for entanglement en- 
tropy by the method of least squares. 



Such distinctively different behaviors of EE for PBC 
and OBC lattices are closely related to VBL formations. 
For example, for the square lattice the leading correction 
to boundary spin correlations and EE is associated with a 
group of elemental valence bond loops mentioned above, 
each connecting two nearest neighbor boundary spins. In 
a lattice with rectangular geometry (OBC), the number 
of such elemental loops is always 1 less than the number 
of boundary sites L y . As a result, the leading correction 
to entropy per boundary spin is proportional to (L y — 
1)/L y , which approaches 1 in the thermodynamic limit. 
Accordingly, we expect that the entanglement entropy 
converges slowly following the function 1/L y , as observed 
in Fig. El 

For a lattice wrapped on a cylinder (PBC), winding- 
loops associated with periodic boundary provide extra 
spin-spin correlations to the system. The EE in a finite 
system is lower than the result in the thermodynamic 
limit, due to the fact that such winding valence bond 
loops are more prominent in a smaller lattice compared 
with those in a larger one. Furthermore, amplitudes as- 
sociated with winding loops decay exponentially fast ac- 
cording to the number of boundary spins L y , which is 
of the same order of winding loops' length. As a result, 
we observe a exponentially converging behavior of the 
entanglement entropy in a PBC lattice. 
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FIG. 3: (Color online) Entanglement spectra of the (left) 
square and (right) hexagonal VBS states with cylindrical ge- 
ometry (PBC). In both cases, L y — 16 and L x — 5, in which 
case results have converged to the two-dimensional limit. The 
ground state energy of He is denoted by Ao- The total spins 
S of the states are marked by different symbols. 



B. Entanglement spectrum 

Let us now consider the ES defined as a set of eigen- 
values of the entanglement Hamiltonian H-& := — hip a- 
In a system with PBC, we are able to study the ES as a 
function of momentum in the y-direction (fc in Fig. [3]). As 
shown in Fig. El the ES obtained from the square (hexag- 
onal) VBS states resembles the spectrum of the spin-1/2 
antiferromagnetic (ferromagnetic) Heisenberg chain. The 
reason for the dependence of the ES on the lattice struc- 
ture is as follows: Although both the square and hexago- 
nal lattices are bipartite, neighboring boundary spins be- 
long to the different sublattices in the case of the square 
VBS state. As a result, the entanglement Hamiltonian is 
reminiscent of the AFM Heisenberg chain. In contrast, 
in the hexagonal lattice, all boundary spins belong to the 
same sublattice, hence the FM spectrum appears. The 
lowest-lying modes in the left panel can be identified as 
the des Cloizeaux-Pearson spectrum in the AFM Heisen- 
berg chain^i. In the right panel, the excitations with 
the total spin S = 7 look like the ordinary spin-wave 
spectrum in the FM Heisenberg chain. In fact, due to 
the translational symmetry in the y direction, the Bloch 
theorem applies and the single-magnon states in the FM 
chain are exact eigenstates of the RDM for the hexagonal 
VBS state. 



C. Nested entanglement entropy 

In order to further establish the correspondence be- 
tween the holographic spin chain and the Hamiltonian 
for the Heisenberg chain, we introduce a measure, which 
we call nested entanglement entropy (NEE), and study its 
scaling properties. One might think that the finite-size 
scaling analysis of the ground state energy Aq is suffi- 



TABLE II: Slopes (0 e av) of modes at k = for L v = 16. 



1.69000 



2.02221 



2.13535 



La, =4 
2.17006 



J^x — 5 
2.18755 



cient to reveal the CFT structure in the holographic spin 
chain. However, there is a subtle point here. The spec- 
trum of the RDM can only tell us that the entanglement 
Hamiltonian can be expressed as He = j3 gH\ 10 \ with 
the Hamiltonian for the holographic chain H\ w \. There 
is no unique way to disentangle a fictitious temperature 
/3 e fi from -ffhoi- Assuming that if hoi is gapless and its 
low-energy dispersion is given by vk, we can estimate 
/3 c ffU from the slope of the modes at k = in the left 
panel of Fig. |3] (see Table HI|). The obtained slope shows 
saturation at about L x = 4, 5, which suggests that the 
fictitious temperature is presumably nonzero (/3 c g < oo) 
in the infinite 2D system. As we will see, however, the 
NEE provides a more clear-cut approach to investigate 
the critical behavior of ifhoi without any assumption. 

Let us now give a precise definition of the NEE. Based 
on the ground state of the entanglement Hamiltonian He, 
we construct the nested RDM for the sub-chain of length 
£ in the holographic spin chain as 



p{£) = Ttt+i 



Mo)(M, 



(12) 



where jV'o) denotes the normalized ground state of He 
and the trace is taken over the remaining sites excluding 
the sub-chain. The NEE is then obtained from the nested 
RDM as 



S(t,L v ) = -Tr li ... tt \p(j()hip(£)], 



(13) 



where the trace is over the sites on the sub-chain. Since 
the low-energy physics of the AFM Heisenberg chain is 
described by c = 1 CFT, the NEE is expected to behave 



s PBC (e,L y ) 



-ln[/(£)] + Sl , 



L v . fn£ 

— - sin — 

7T \L V 



(14) 
(15) 



for a lattice with PBC, where c is the central charge, and 
si is a non-universal constant. In Fig. U (a), we show 
NEE obtained from the square VBS state with L x = 



5 and L y = 16. The CFT prediction, Eq. flH}, well 
explains the spatial profile of the NEE. The fit yields 
c = 1.01(7) which is reasonably close to c = 1. 

The NEE for I = L y /2 (half-NEE) is simplified to be 
S PBC (L v /2,L y ) = (c/3)ln(L y ) + const. From the data 
with the finite-size scaling form of the NEE, we can also 
extract the central charge c as summarized in Tabic Hill 
For the case L x = 1, c is remarkably close to unity. This 
further supports our conjecture: a correspondence be- 
tween the entanglement Hamiltonian of the 2D AKLT 
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FIG. 4: (Color online) Nested entanglement entropy S(£, L y ) 
as a function of the sub-chain length £ for L x — 5 and L y — 16. 
(a) and (b) show results obtained for square VBS states with 
PBC and OBC, respectively. Fits to the CFT predictions, 
Eqs. (|14[) and (|16[) . are indicated by open circles. The lines 
are guides to the eye. 



TABLE III: The obtained central charge from half-NEE. 



L m = l 



h x — 2i 



L x — o 



L x =4 



L*x — o 



1.007(4) 1.042(4) 1.055(4) 1.056(2) 1.059(2) 



model and the physical Hamiltonian of the ID Heisen- 
berg chain. A slight modification of c is observed as the 
bulk width increases from L x = 1 to L x = 5. 

In the case of the rectangular geometry (OBC), a stag- 
gered pattern of NEE is expected from studies of the 
standard EE in open spin chain o 26 ' 27 : 



S°^(£,L y ) = -h%(2£+l)] + a 



7TCl 



i-iy 



M = 



G 

2(i, 



2v 5 (2£+l) 1 /2i<" 



1) 



sin 



nf 



2{L y 



(16) 



In Fig. 0] (b), we show the NEE for the square VBS state 
on the rectangle with L x = 5 and L y = 16. We take the 
central charge to be c = 1 and the Tomonaga-Luttingcr 
parameter K = 1, and tune a and C\jv as fitting param- 
eters. The fit is reasonably good except for the bound- 
aries. This deviation might be due to logarithmic cor- 
rections which appear in the SU(2) symmetric AFM spin 
chains 2 ^ but are not taken into account in Eq. (JTHJ). 

To provide further evidence for the holographic spin 
chain, we now consider the hybrid system comprised of 
squares and hexagons, shown in Fig.[5{a). It is expected 
that the entanglement Hamiltonian is described by the 
FM-AFM alternating Heisenberg model. Figure [Sfb) 
shows the ES for this system. It is clearly seen that 
there is a gap between the ground state and the contin- 
uum. The gap rapidly saturates with increasing L y . The 
behavior of the ES is totally consistent with the energy 
spectrum of the alternating Heisenberg chain whose in- 
teractions are denoted J and J'. 28 Comparing the bound- 
state spectrum in Fig. [SJb) with that of the triplet-pair 
excitation in the alternating Heisenberg chain^ we esti- 
mate the ratio of the FM exchange divided by the AFM 




the reduced density matrix of a subsystem can be inter- 
preted as a thermal density matrix of the holographic 
spin chain, whose spectrum resembles that of the spin- 
1/2 Hcisenberg chain. To elucidate this relationship, we 
have introduced the concept, nested entanglement en- 
tropy (NEE), which allows us to clarify this correspon- 
dence in a quantitative way without information on the 
fictitious temperatures. The finite-size scaling analysis of 
the NEE revealed that the low-energy physics of the holo- 
graphic chain associated with the square lattice VBS is 
well described by c = 1 conformal field theory. The NEE 
was also applied to the hybrid VBS state, the entangle- 
ment spectrum of which is gapped, and the holographic 
chain was found to be the alternating Heisenberg chain. 
Acknowledgment 



FIG. 5: (Color online) (a) (Left) VBS state on a hybrid 
lattice. The red (dotted) line indicates the reflection axis. 
(Right) Corresponding holographic chain. The thick and 
dashed lines indicate AFM and FM bonds, respectively, (b) 
Entanglement spectrum of the VBS state in (a) with L y = 16. 
The ground state energy is denoted by Ao. The dashed curve 
indicates the spectrum of triplet-pair bound states, (c) Nested 
entanglement entropy S(£, L y ) as a function of £ for L y = 16. 



one as J'/ J ~ 0.5. Note that one can further manipulate 
J' I J of the holographic spin chain by tuning the lattice 
structure and the number of valence bonds on each edge 
the lattice. The NEE associated with the ground state 
of H E is shown in Fig. EJc). The NEE is 2 log 2 when 
two AFM bonds are cut, whereas it becomes about log 2 
when both FM and AFM bonds are cut. The obtained 
result is in good agreement with the standard EE in the 
alternating Heisenberg chain studied in Ref. [29|. 



IV. CONCLUSION 

We have studied both the entanglement entropy and 
spectrum associated with the VBS ground state of the 
AKLT model on various 2D lattices. It was shown that 
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Appendix A: Analytical method for L x = 1 

In this appendix, we show how to explicitly calculate 
the reduced density matrix for L x — 1 ladders with OBC 
and PBC based on an extension of the method described 
in Ref. [22|, in which the authors analytically obtained 
entanglement entropies of the VBS states with OBC us- 
ing a transfer matrix technique. According to Ref. |22| . 
the overlap matrix M for vertical ladder, i.e. L x = I is 
given by 



J 



M 




L y - 1 (OBC) 
(PBC) 



(Al) 



where zi denotes the coordination number at the i- 
th site, Qi is the unit vector defined by f2^ := 
(sin 6^ cos </> , , sin j sin </>j , cos 0, ) , and Bj = (a^a^ tr|) 

is the Pauli matrix vector. For PBC, (II +i = &i. 
The above expression for M can be derived using the 
Schwinger boson representation and coherent state repre- 



sentation. It is convenient to introduce Z-matrix defined 



by 



Z(L v ;{a i} ) :-/fn^jn(l + A*-S! 



Jl(i-n k -Ci k+1 ), 



(A2) 



fc=i 



Since the only difference between M and Z is the overall 
factor, the spectrum of pa is equivalent to that of the 
following matrix: 



Pa 



z{l v ;{^}) 2 

Tr[Z(L y ;{a t }n 



(A3) 



where Tr is taken over tr-spin spaces. Then in order to 
calculate entanglement properties of the VBS state on a 
vertical ladder, it is necessary to obtain the explicit form 
of Z-matrix. There is a useful relation shown in lemma 
3.3 in Ref. fil: 



-7-^(<?i ■ &i){&i ■ crj) = q&i 



(A4) 



where q = 1/3. Using the above relation we obtain the 
following relation: 



Z(L y ;{ai}) 



X (o»i • Vi 2 ) ■ ■ • (<?i2n- 



(-«)' 



(A5) 



where lj = mod(z2j — iij-x, L y ) and 1 < i\ < ii < ■ ■ ■ < 
*2n— l < Ly. The summation is taken over admissible 

combinations of {ij}j = i 2n for < n < N/2. It is 

useful to express Eq. (|A5|) in terms of graphical repre- 
sentations. 

Figure [5] shows the relationship between the origi- 
nal VBS state and the one-dimensional holographic spin 
chain. The circles in the right panel in Fig.[S]indicate the 
sites in the spin chain, on which the Pauli vectors {ffi} 




FIG. 6: (Color online) Relationship between the original 
VBS state on a square-lattice vertical ladder (left) and a holo- 
graphic spin chain (right). The circles in the right panel in- 
dicate the sites in the spin chain, on which the Pauli vectors 
{ffi} in Eq. (KE\i act. 



in Eq. (|A5[) act. Z-matrix is represented by sum of ad- 
missible configurations each of which has a corresponding 
graph (see Fig. [7]). Here the number of admissible con- 
figurations is 2 Ly . 

Let us consider the L y = 4 case for illustrative pur- 
poses. Figure UJa) shows all admissible configurations 
for OBC using a bit representation in which the thin and 
thick lines denote and 1, respectively. The end points 
of thick lines in Fig. [7] indicate {ij}. For example, the 
contributions of the states 3 and 5 to the Z-matrix are 
expressed as 



q 1 (-q) 2 (<7i -a 3 ), 
<?{-q) l {-<i) 1 ($1 ■0 ; 2)(<?3 -Si) 



(A6) 
(A7) 



respectively. Note that n in Eq. (|A5|) corresponds to the 
number of clusters connected by thick lines contiguously 
in each graphical representation. Summing all contribu- 
tions, we arrive at the following expression for the Z- 
matrix: 

Z OBC (4; {«?,}) =I-q 2 [(*i • ff 2 ) + {32 ■ ff 3 ) + (<?3 ■ <? 4 )] 

+q 3 [((Ji • <T 3 ) + (<7 2 • 5i)] 

W [(ffl ■ <? 2 ) ((7 3 • (7 4 ) - (ffi -ffi)], (A8) 

where I is a 2 Ly -dimensional identity matrix. 

Next, we examine the case of PBC where the number 
of admissible configurations is 2 Ly (— 2 Ly ). We have to 
consider the states depicted in Fig. [7] (b) in addition to 
the states in Fig. [JJ(a). For example, contribution of the 
state 11 to the Z-matrix is given as 



q l {-qf{a A -aj,). 



(A9) 



Including additional contributions from the states from 
8 to 15 in Fig. [7J(b), the Z-matrix for PBC is given as 

Z PBC (4; {<?,}) =Z OBC (4;{(7j)- g 2 (oWi) 
+q 3 [(ffi ■ ff 2 ) + {ffa ■ ffi) + I] 
+q 4 [((7 4 • ffi) (ct 2 • era) - (<? 4 • ffa) 
- ((7 3 • ff 2 ) - {ffa ■ ffi)] ■ (A10) 

In principle, we can obtain an explicit form of the Z- 
matrix for larger L y using the method based on the 
graphical representations. 

Next we show how to calculate the Z-matrices for 
hexagonal-lattice vertical ladders (i.e. L x — 1). The 
relationship between the VBS on the hexagonal-lattice 
vertical ladder and the holographic spin chain is shown 
in Fig. [8) The only difference between the square and 
hexagonal vertical ladders is the definition of the length 
between sites on the holographic spin chain. We can then 
calculate the Z-matrix for the hexagonal-lattice vertical 
ladder as follows: 



Z{L y -{ff t }) = E 1 n (-l) 2h 

{il-/-- ,t2n} 

•l) ' ' ' { a i2n-l 



{-if 



X [<Ji 



(AH) 



3 4 
(a) 




FIG. 7: (Color online) Graphical representation of all admissible graphs. When we study the OBC case, we have to consider 
only the configurations depicted in (a), whereas we have to consider the configuration depicted in (b) in addition to that 
depicted in (a), we study the PBC case. 




FIG. 8: (Color online) Relationship between the original 
VBS state on a hexagonal-lattice vertical ladder (left) and a 
holographic spin chain (right). The circles in the right panel 
indicate the sites in the spin chain, on which the Pauli vectors 
{ffi} in Eq. (KE\i act. 



where lj = modfaj — i2j-i,L y ) and 1 < i\ < 12 < 
■ ■ ■ < i2n-i < L y . The difference between Eq. (|A5[) and 
Eq. (jAlip is an exponent of (—<?). 

After we obtain an explicit form of Z-matrix, we can 
study entanglement properties. Indeed, we have com- 
pared the results based on this method with those ob- 
tained by Monte Carlo method and have confirmed that 
they are completely consistent, which ensures the validity 
of our Monte Carlo approach. 



Appendix B: Monte Carlo Method for Calculating 
Overlap Matrix 

In this appendix, we explain our proposed Monte Carlo 
method used for obtaining the overlap matrix M. As we 
mentioned in Sec. II, the VBS state can be regarded as a 
singlet-covering state, expressed by Eq. @. In the case 
of bipartite lattices (square and hexagonal), we can avoid 
negative signs by a local gauge transformation: 



„(i)t 



kt 



(ieA) 



„(2)t 



,t 



(ieA) 
(i e B) 



(Bl) 



where A and B denote the two sublattices. 
state can then be represented as 



The VBS 



l*> = 



IK 



c (l)t c (Dt 



„(2)f „(2)t 



)|vac), 



(B2) 



where (i, j) runs over all connected neighboring sites in 
the graph. If we introduce r = 1, 2 as the quantum num- 
ber of SU(2) boson operators, the VBS state can be writ- 
ten as 



/ 



I*) 



{i-jk} 



E 114 



{T 3 k)\ (r jk )\ 



\ 



\ U,k) 

\j,fcefi 



I vac). 



(B3) 



/ 



where Tjk denotes the boson between j-th and fc-th site 
on the bulk of the subsystem tt. As a result, we can re- 
produce VBS state through stochastic sampling of boson 
configurations {rjk} of the system. 

The element of overlap matrix M , is given by 

M W},{-.} = (W}IW>. ( B4 ) 

where |{ci}) is the VBS state with fixed boundary {<Ji}: 

\ 



{r ]k } Viean 

= E li^'i 7 "^})' 
i T Jk} 



n 

(l,k) 



„( T ifc)t^,(Tjfc)t 



/ 



I vac) 



(B5) 



where <9f2 denotes the boundary of $7. The variables 
{<7i},i € dQ specifies the unpaired spins on the bound- 
ary. We can further simplify the problem by considering 
(IB5P in the occupation number basis: 



{ T jk} 



(B6) 



where n a = ^2 ke g .6(Tjk,<x) (the summation is taken 
over all nearest neighbor sites of j) is the number of 
bosons of type a at the site j. The degeneracy of oc- 
cupation number state can be easily derived as: 



w({rij }) = ]_ j v a 'a 



(l), n (2), 



(B7) 



As a result, we obtain the element of the overlap matrix 



M. 



{<},{°i} 



= E W({n^})5({n ( p k, },{n^}), 



{T' ]k },{T 3 k} 



(B8) 



with weight 



(r jk ) x , _ 



W({nf 3k '})=w 2 (nY 3k ') 



(Tjk)\ _ 



ik 



Ml (2),\ 

rij '\n) '\) 



(B9) 



In the end, by generating boson states £ = 

({r' fe }, {Tjfe}) stochastically while keeping the constraint 

( T '\ ( T \ 
n- = n- for all j and r, we can numerically calculate 

the element of the overlap matrix Ms a >\ t ai \ as accumula- 
tion number of boson configurations with corresponding 
boundary condition, through Monte Carlo importance 
sampling. It is convenient to introduce graphical repre- 
sentation for each boson configurations |{ci}){ryfc}), as 
shown in Fig. ® The red (dark) and blue (light) lines 
indicate bosons with quantum number Tjk — 1 and 2, 
respectively. Colored but unconnected dots on the edge 
correspond to boundary bosons {<Ji\- Notice that we 
show two layers in the graphical representation, as they 
stand for the sub-system state I^I" 4 ]) and its conjugate 
(\JjMl, respectively. Due to the fact that the full system 
is reflection symmetric, one may also consider two lay- 
ers as two subsystems A and B, with sites (big circles) 
and bonds (lines) on top of each other being symmetric 
counterparts. 



Oi (Ti\ O 







m 5 



GfOG 



FIG. 9: (Color online) Illustration of boson configurations 
(states) and two MC update schemes: (left) Inter-layer flip, 
(right) intra-layer flip. Red (dark) and blue (light) lines stand 
for bosons Tjk = 1,2, respectively. Upper and lower layers 
stand for a state in sub-systems A (|\I/' ')) and its conjugate 
((*&' '|), respectively. Valence bonds being flipped are indi- 
cated by a yellow (shadowed) plaquette. The right edge of 
each panel corresponds to the boundary between the subsys- 
tems A and B, and the unconnected dots represent the a 
variables. 



Using this graphical representation, we can easily il- 
lustrate two MC update schemes we use to sample all 
possible boson configurations, while automatically satis- 



fying the constrain: n- 



(r') 



» 



for all j and r. 



(i) Inter-layer flip (Left panel of Fig. |H|): In the boson 
state E, we choose arbitrarily a valence bond connect- 
ing sites j and k such that Tjk = r'- k = fj, (colored as 
red/blue) is satisfied. A trial state St is given as Tjk — 
tL = 3 — fi (colored as blue/red) with all other bosons 
unchanged. The update is accepted/rejected according 
to the ratio of probabilities W(Et)/W(£) which satisfies 
detailed-balance condition. In Fig. |H1 a yellow loop is 
drawn showing bosons being flipped (Tjk = H — > 3 — /i). 

(ii) Intra-layer flip (Right panel of Fig. [9]): In order to 
apply intra-layer flip without violating the constraints, 
we attempt to construct a "valence bond loop" (VBL), or 
a sequence of connected singlets, with alternating colors 
(Tjk = 1h2). Flipping a VBL is equivalent to operating 



n 



(a„)t (a,.)t (A.) (A.) 



« K =/w = r , { ; : odd) j (bio) 

I v = 3 — fi (k : even) / 



on the VBS state (except for the weight change), where 
the loop is formed by / sites {i K } in the order of k = 
1,2,...,/, and fi,v are the alternating colors along the 
loop, (from red to blue and vice versa). In the right col- 
umn of Fig. El we show a possible loop that consists of 4 
sites and 4 valence bonds, denoted by yellow region. In 
our MC simulation scheme, we construct a loop by using 
a random walker. The walker starts from a randomly 
chosen site, as well as a randomly chosen direction. Af- 
ter such a movement, we choose next direction so that 
the color of next valence bond alternates. Since turn- 
ing back to the previous site is in the list of candidates, 
it is always possible to choose such a direction. When 
the walker comes back to the initial position, the loop 
is formed. Therefore, loops can be of arbitary length, 
as long as color is alternating along the loop. Mean- 
while, the possibility (amplitude) decrease significantly 
for longer loops. Nevertheless, obviously, flipping all the 
bonds in a valence bond loop does not affect the occupa- 
tion number n- a at any site on the loop. 

Following above two updating rules, we can obtain effi- 
ciently the overlap matrix with high accuracy, while sat- 
isfying the ergodicity condition. This MC sampling ap- 
proach can be applied to arbitrarily large systems. Num- 
ber of sites along the bulk dimension has limited effect 
on the cost of calculations. The computational limit is 
mainly related to the diagonalization of the overlap ma- 
trix, the dimension of which grows exponentially as the 
number of boundary sites increases. 
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